Integrative Analysis

Load Libraries & ArchR object

library(ArchR)
library(rhdf5)
library(reticulate)
library(tidyverse)
#library(caret)
h5disableFileLocking()

# load ArchR project
proj <- loadArchRProject(path = "Save_Vignette_object_4/")

# for R to find macs2 we will use a conda environment where the latest version 
# of macs2 is installed
#conda_list()[[1]][2] %>% use_condaenv(required = TRUE)

Create Low-Overlapping Aggregates of Cells

If we compute correlation on sparse features, this can lead to noise in the data. Therefore, first low-overlapping aggregates of single cells are created (see Cicero). Aggregates with >80% overlap with an other aggregate are removed.

• sampling with replacement, k = 100 (default) groups/cell
• aggregate binary accessibility profiles for cells in each group

Co-accessibility

Correlation in accessiblity between two peaks across many single cells. For example if peak A is accessible in a single cell is peak B also accessible? With this analysis we might learn if a certain enhancer peak is often co-accessible with a certain promoter peak in a particular cell type.

Keep in mind that cell-type specific peaks would be identified as co-accessible within this cell tpye. This strong correlation does not mean that there is an underlying co-accessibility!

Query hits and subject hits are the indices of the two peaks that are correlated. The correlation column contains the numeric correlation value of the accessibility between the two peaks.

# compute peak co-accessibility information and store it in the ArchR project
proj <- addCoAccessibility(proj, 
                           reducedDims = "IterativeLSI",
                           k = 100,
                           corCutOff = 0.75,
                           # max. 80 % overlap between current group and all previous groups to permit
                          # current grup to be added to the group list during k-nearest neighbor calculation
                          overlapCutoff = 80,
                          dimsToUse = 1:30,
                           #number of k-nearest neighbor groupings to test for passing the supplied overlapCutoff
                          knnIteration = 500 )

# to retrieve the co-accessibility information we use getCoaccessibility()
# return dataframe if returnLoops = FALSE
coacc <- getCoAccessibility(
    ArchRProj = proj,
    # if the dimension has a correaltion to the sequencing depth 
    # that is greater than the cutoff it will be excluded
    corCutOff = 0.5,
    resolution = 1,
    returnLoops = FALSE # if TRUE it returns a GRangesObject
)

coacc %>% head %>% knitr::kable()

queryHits	subjectHits	seqnames	correlation	Variability1	Variability2	TStat	Pval	FDR	VarQuantile1	VarQuantile2
6	9	chr1	0.5948088	0.0041591	0.0302828	16.51231	0	0	0.4991376	0.9302221
9	6	chr1	0.5948088	0.0302828	0.0041591	16.51231	0	0	0.9302221	0.4991376
9	31	chr1	0.6909559	0.0302828	0.0141892	21.32989	0	0	0.9302221	0.8024554
18	19	chr1	0.5128650	0.0259368	0.0052909	13.33193	0	0	0.9060243	0.5682542
19	18	chr1	0.5128650	0.0052909	0.0259368	13.33193	0	0	0.5682542	0.9060243
21	49	chr1	0.5345478	0.0042255	0.0335266	14.11476	0	0	0.5038296	0.9449599

Metadata:

metadata(coacc)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
##             seqnames              ranges strand
##                <Rle>           <IRanges>  <Rle>
##        Mono     chr1       752506-753006      *
##           B     chr1       762688-763188      *
##         GMP     chr1       773653-774153      *
##           B     chr1       801007-801507      *
##           B     chr1       805041-805541      *
##         ...      ...                 ...    ...
##        Mono     chrX 154840753-154841253      *
##          NK     chrX 154842390-154842890      *
##          NK     chrX 154862033-154862533      *
##   Erythroid     chrX 154912392-154912892      *
##         GMP     chrX 154997007-154997507      *
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
dim(coacc)
## [1] 111132     11

We can also return the co-accessibility information as a loop track. The resolution parameter sets the base-pair resolution of the loops. Resolution = 1 creates loops that connect the center of each peak.

Where is the information contained which peaks are correlated?

coacc <- getCoAccessibility(ArchRProj = proj,
                            corCutOff = 0.5,
                            resolution = 1,
                            returnLoops = TRUE)


coacc [[1]]
## GRanges object with 55566 ranges and 9 metadata columns:
##           seqnames              ranges strand | correlation Variability1
##              <Rle>           <IRanges>  <Rle> |   <numeric>    <numeric>
##       [1]     chr1       845703-856628      * |    0.594809   0.00415906
##       [2]     chr1       856628-940499      * |    0.690956   0.01418923
##       [3]     chr1       894695-895356      * |    0.512865   0.00529094
##       [4]     chr1       901579-999702      * |    0.534548   0.00422549
##       [5]     chr1      948868-1004761      * |    0.506462   0.02580794
##       ...      ...                 ...    ... .         ...          ...
##   [55562]     chrX 153583564-153665944      * |    0.537871  0.000662602
##   [55563]     chrX 153584852-153597227      * |    0.566064  0.001244477
##   [55564]     chrX 153597227-153637613      * |    0.512202  0.063475852
##   [55565]     chrX 153686270-153744595      * |    0.513243  0.012864089
##   [55566]     chrX 153947827-154027599      * |    0.524475  0.013135018
##           Variability2     TStat        Pval         FDR VarQuantile1
##              <numeric> <numeric>   <numeric>   <numeric>    <numeric>
##       [1]    0.0302828   16.5123 3.62193e-49 2.50640e-47     0.499138
##       [2]    0.0302828   21.3299 3.41000e-72 6.81189e-70     0.802455
##       [3]    0.0259368   13.3319 6.84909e-35 2.23707e-33     0.568254
##       [4]    0.0335266   14.1148 2.71023e-38 1.06988e-36     0.503830
##       [5]    0.0189477   13.1076 6.23013e-34 1.92659e-32     0.905221
##       ...          ...       ...         ...         ...          ...
##   [55562]   0.00116998   14.2381 7.75355e-39 3.14586e-37    0.0623195
##   [55563]   0.06347585   15.3237 1.05928e-43 5.56776e-42    0.1644870
##   [55564]   0.02888357   13.3086 8.62737e-35 2.80265e-33    0.9912758
##   [55565]   0.02225960   13.3453 6.00370e-35 1.96741e-33    0.7838576
##   [55566]   0.00519124   13.7465 1.10713e-36 3.99906e-35    0.7879772
##           VarQuantile2     value
##              <numeric> <numeric>
##       [1]     0.930222  0.594809
##       [2]     0.930222  0.690956
##       [3]     0.906024  0.512865
##       [4]     0.944960  0.534548
##       [5]     0.853012  0.506462
##       ...          ...       ...
##   [55562]     0.151581  0.537871
##   [55563]     0.991276  0.566064
##   [55564]     0.923357  0.512202
##   [55565]     0.880343  0.513243
##   [55566]     0.562910  0.524475
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

If we decrease the resolution to 1000 this could help with overlplotting of co-accessibility interactions. The GRanges object will contain fewer total entries, 51,82 instaed of 54,002.

cA <- getCoAccessibility(
    ArchRProj = proj,
    corCutOff = 0.5,
    resolution = 1000,
    returnLoops = TRUE
)

cA[[1]]
## GRanges object with 53116 ranges and 9 metadata columns:
##           seqnames              ranges strand | correlation Variability1
##              <Rle>           <IRanges>  <Rle> |   <numeric>    <numeric>
##       [1]     chr1       845500-856500      * |    0.594809   0.00415906
##       [2]     chr1       856500-940500      * |    0.690956   0.01418923
##       [3]     chr1       894500-895500      * |    0.512865   0.00529094
##       [4]     chr1       901500-999500      * |    0.534548   0.00422549
##       [5]     chr1      948500-1004500      * |    0.506462   0.02580794
##       ...      ...                 ...    ... .         ...          ...
##   [53112]     chrX 153583500-153665500      * |    0.537871  0.000662602
##   [53113]     chrX 153584500-153597500      * |    0.566064  0.001244477
##   [53114]     chrX 153597500-153637500      * |    0.512202  0.063475852
##   [53115]     chrX 153686500-153744500      * |    0.513243  0.012864089
##   [53116]     chrX 153947500-154027500      * |    0.524475  0.013135018
##           Variability2     TStat        Pval         FDR VarQuantile1
##              <numeric> <numeric>   <numeric>   <numeric>    <numeric>
##       [1]    0.0302828   16.5123 3.62193e-49 2.50640e-47     0.499138
##       [2]    0.0302828   21.3299 3.41000e-72 6.81189e-70     0.802455
##       [3]    0.0259368   13.3319 6.84909e-35 2.23707e-33     0.568254
##       [4]    0.0335266   14.1148 2.71023e-38 1.06988e-36     0.503830
##       [5]    0.0189477   13.1076 6.23013e-34 1.92659e-32     0.905221
##       ...          ...       ...         ...         ...          ...
##   [53112]   0.00116998   14.2381 7.75355e-39 3.14586e-37    0.0623195
##   [53113]   0.06347585   15.3237 1.05928e-43 5.56776e-42    0.1644870
##   [53114]   0.02888357   13.3086 8.62737e-35 2.80265e-33    0.9912758
##   [53115]   0.02225960   13.3453 6.00370e-35 1.96741e-33    0.7838576
##   [53116]   0.00519124   13.7465 1.10713e-36 3.99906e-35    0.7879772
##           VarQuantile2     value
##              <numeric> <numeric>
##       [1]     0.930222  0.594809
##       [2]     0.930222  0.690956
##       [3]     0.906024  0.512865
##       [4]     0.944960  0.534548
##       [5]     0.853012  0.506462
##       ...          ...       ...
##   [53112]     0.151581  0.537871
##   [53113]     0.991276  0.566064
##   [53114]     0.923357  0.512202
##   [53115]     0.880343  0.513243
##   [53116]     0.562910  0.524475
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

Plot Browser Tracks of Co-accessibility

We can add the co-accessibility information as loop tracks to our browser track. For this we can use the loops parameter.

# marker genes we are interested in 
markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

# plot
p <- plotBrowserTrack(
    ArchRProj = proj, 
    groupBy = "Clusters2", 
    geneSymbol = markerGenes, 
    upstream = 50000,
    downstream = 50000,
    loops = getCoAccessibility(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
##       seqnames              ranges strand |     gene_id      symbol
##          <Rle>           <IRanges>  <Rle> | <character> <character>
##   [1]     chr1 208059883-208084683      - |         947        CD34
##   [2]     chrX   48644982-48652717      + |        2623       GATA1
##   [3]     chr9   36838531-37034476      - |        5079        PAX5
##   [4]    chr11   60223282-60238225      + |         931       MS4A1
##   [5]     chr5 140011313-140013286      - |         929        CD14
##   [6]    chr11 118209789-118213459      - |         915        CD3D
##   [7]     chr2   87011728-87035519      - |         925        CD8A
##   [8]    chr17   45810610-45823485      + |       30009       TBX21
##   [9]     chr5   35856977-35879705      + |        3575        IL7R
##   -------
##   seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

grid::grid.newpage()
grid::grid.draw(p$CD14)

Peak2Gene Linkage

What is the difference between co-accessibility and peak2gene linkage?

Co-accessibility uses only scATAC-seq data looking for correlations in accesssiblity.

Peak2Gene linkage in comparison also uses integrated scRNA-seq data to look for correlations between peak accessibility and gene expression.

Both approaces are solutions to a similar problem, but peak2gene linkage is assumed to provide more biologically relevant regulatory links.

All potential peak-to-gene linkages are identified and significant links (R > 0.45 and FDR < 0.1) are kept (Higher correlation and variance suggest that a link is not random noise.)

proj <- addPeak2GeneLinks(
  ArchRProj = proj,
  reducedDims = "IterativeLSI",
  useMatrix = "GeneIntegrationMatrix",
  dimsToUse = 1:30,
  k = 100, # number of k-nearest neighbors to use for creating single cell goups for correlation analysis
  maxDist = 250000, # maximum allowable distance in basepairs between two peaks to consider for co-accessiblity
  scaleTo = 10^4 # total insertion coutns from a group of cell is summed across all peaks and normalized to toal depth provided by scaleTo
)

p2g <- getPeak2GeneLinks(
    ArchRProj = proj,
    corCutOff = 0.45,
    resolution = 1,
    returnLoops = FALSE,
    varCutOffATAC = .25, # minimum variance quantile of ATAC peak accessibility
    varCutOffRNA = .25, # minimum variance quantile of the RNA gene expression when selecting links
    FDRCutOff = 1e-04 # maximum numeric peak2gene FDR to return
)

as_data_frame(p2g) %>% arrange(desc(Correlation)) %>% 
  head %>% knitr::kable(caption = "Peak-to-gene linkage")
## Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
## Please use `as_tibble()` instead.
## The signature and semantics have changed, see `?as_tibble`.

Peak-to-gene linkage

idxATAC	idxRNA	Correlation	FDR	VarQATAC	VarQRNA
76211	2257	0.9776989	0	0.8779933	0.6813612
42485	14306	0.9621724	0	0.9966034	0.9910757
23457	10826	0.9599931	0	0.9372381	0.5969572
81836	3162	0.9599560	0	0.9999645	0.9820440
30183	12832	0.9559511	0	0.2635277	0.5640019
7959	1140	0.9556089	0	0.9458961	0.3378851

In the metadata the peakSet with ranges, strand and chromosome are saved, as well as the gene set with ranges, strand, chromosome and gene name.

metadata(p2g)
## $peakSet
## GRanges object with 141025 ranges and 0 metadata columns:
##            seqnames              ranges strand
##               <Rle>           <IRanges>  <Rle>
##        [1]     chr1       752506-753006      *
##        [2]     chr1       762688-763188      *
##        [3]     chr1       773653-774153      *
##        [4]     chr1       801007-801507      *
##        [5]     chr1       805041-805541      *
##        ...      ...                 ...    ...
##   [141021]     chrX 154840753-154841253      *
##   [141022]     chrX 154842390-154842890      *
##   [141023]     chrX 154862033-154862533      *
##   [141024]     chrX 154912392-154912892      *
##   [141025]     chrX 154997007-154997507      *
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
## 
## $geneSet
## GRanges object with 18601 ranges and 2 metadata columns:
##           seqnames    ranges strand |      name     idx
##              <Rle> <IRanges>  <Rle> |   <array> <array>
##       [1]     chr1     69091      * |     OR4F5       1
##       [2]     chr1    762902      * | LINC00115       2
##       [3]     chr1    812182      * |    FAM41C       3
##       [4]     chr1    860530      * |    SAMD11       4
##       [5]     chr1    894679      * |     NOC2L       5
##       ...      ...       ...    ... .       ...     ...
##   [18597]     chrX 154299695      * |     BRCC3     708
##   [18598]     chrX 154444701      * |      VBP1     709
##   [18599]     chrX 154493852      * |    RAB39B     710
##   [18600]     chrX 154563986      * |     CLIC2     711
##   [18601]     chrX 154842622      * |     TMLHE     712
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths
## 
## $seATAC
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seATAC-Group-KNN.rds"
## 
## $seRNA
## [1] "/omics/groups/OE0533/internal/katharina/scDoRI/ArchR/Save_Vignette_object_4/Peak2GeneLinks/seRNA-Group-KNN.rds"

atac_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seATAC-Group-KNN.rds")
rna_knn <- readRDS("Save_Vignette_object_5/Peak2GeneLinks/seRNA-Group-KNN.rds")
atac_knn
## class: RangedSummarizedExperiment 
## dim: 141025 491 
## metadata(1): KNNList
## assays(2): ATAC RawATAC
## rownames(141025): f1 f2 ... f141024 f141025
## rowData names(13): score replicateScoreQuantile ... idx N
## colnames: NULL
## colData names(0):

dim(assays(atac_knn)[[1]])
## [1] 141025    491
#assays(atac_knn)[[1]] %>% head

#assays(rna_knn)[[1]] %>% head

We can also return loops:

p2g_loop <- getPeak2GeneLinks(
    ArchRProj = proj,
    corCutOff = 0.45,
    resolution = 1000, # decrease resolution for plotting
    returnLoops = TRUE
)

p2g_loop[[1]] %>% head
## GRanges object with 6 ranges and 2 metadata columns:
##       seqnames         ranges strand |     value         FDR
##          <Rle>      <IRanges>  <Rle> | <numeric>   <numeric>
##   [1]     chr1  762500-812500      * |  0.450391 7.66477e-25
##   [2]     chr1  762500-895500      * |  0.460307 4.97788e-26
##   [3]     chr1  762500-948500      * |  0.515106 2.34837e-33
##   [4]     chr1  812500-894500      * |  0.451671 5.41109e-25
##   [5]     chr1  812500-895500      * |  0.496871 9.11031e-31
##   [6]     chr1 812500-1004500      * |  0.460707 4.44943e-26
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

markerGenes  <- c(
    "CD34", #Early Progenitor
    "GATA1", #Erythroid
    "PAX5", "MS4A1", #B-Cell Trajectory
    "CD14", #Monocytes
    "CD3D", "CD8A", "TBX21", "IL7R" #TCells
  )

p <- plotBrowserTrack(
    ArchRProj = proj, 
    groupBy = "Clusters2", 
    geneSymbol = markerGenes, 
    upstream = 50000,
    downstream = 50000,
    loops = getPeak2GeneLinks(proj)
)
## GRanges object with 9 ranges and 2 metadata columns:
##       seqnames              ranges strand |     gene_id      symbol
##          <Rle>           <IRanges>  <Rle> | <character> <character>
##   [1]     chr1 208059883-208084683      - |         947        CD34
##   [2]     chrX   48644982-48652717      + |        2623       GATA1
##   [3]     chr9   36838531-37034476      - |        5079        PAX5
##   [4]    chr11   60223282-60238225      + |         931       MS4A1
##   [5]     chr5 140011313-140013286      - |         929        CD14
##   [6]    chr11 118209789-118213459      - |         915        CD3D
##   [7]     chr2   87011728-87035519      - |         925        CD8A
##   [8]    chr17   45810610-45823485      + |       30009       TBX21
##   [9]     chr5   35856977-35879705      + |        3575        IL7R
##   -------
##   seqinfo: 24 sequences from hg19 genome
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

grid::grid.newpage()
grid::grid.draw(p$CD14)

Heatmap of peak2gene linkage

The rows are clustered using k-means clustering based on value passed to parameter k (default = 25)

plotPeak2GeneHeatmap(proj, groupBy = "Clusters2")
## Warning: did not converge in 10 iterations